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I. INTRODUCTION 



Neutrino experiments are closing in on neutrinos in the ultra-high energy (UHE) regime, where a diffuse 
neutrino flux, first predicted by Berezinsky and Zatsepin [1, 2], is expected to result from interactions 
between UHE cosmic rays and cosmic microwave background photons through what is known as the Greisen, 
Zatsepin and Kuzmin (GZK) process [3, 4]. Neutrinos in this energy regime probe higher center-of-mass 
(COM) energies than those accessible by human-made accelerators through their interactions in the earth. 
For example, the COM energy of a fO 9 GeV neutrino incident on a nucleon at rest is 45 TeV. Thus, a 
measurement of neutrino-nucleon {vN) cross sections in the UHE regime could be sensitive to either new 
physics scenarios such as extra space-time dimensions or unexpected behavior of parton distribution functions 
(PDFs) at Bjorken-x smaller than that accessible by current experiments (from here on, any v refers to both 
neutrino and anti- neutrino unless otherwise stated) [5]. However, before the significance of any vN cross 
section measurement can be assessed, the uncertainties on the Standard Model (SM) expectation must be 
quantified based on the diverse body of current experimental constraints. 

The paper is composed of two parts. In Section II, we perform a new calculation of vN cross sections 
and their associated PDF uncertainties for neutrino energies E v > 10 4 GeV using the MSTW 2008 PDF 
set. In Section II A, we review the expressions for the vN cross sections in terms of the quark PDFs. In 
Section II B we present the results of our cross section calculations, associated uncertainties and their energy 
dependent parametrizations. Next, in Section II C we discuss the differential cross sections, and parametrize 
the inelasticity distributions in an energy dependent way. We also show a few select distributions in a;. In 
Section II D, we calculate the correlations between the uncertainties across energies. 

In Section III, we propose to constrain the UHE vN cross section in future sub-terranean neutrino experi- 
ments using a technique that is independent of the incident flux spectrum of neutrinos and present projected 
constraints on models with enhanced cross sections due to extra space-time dimensions. 



II. NEUTRINO-NUCLEON CROSS SECTION 



A. Methodology 

The vN cross section for charged current (CC) interactions on an isoscalar target is given by 1 : 
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with quark and antiquark densities given by q = (d + u)/2 + s + b and q = (d + u)/2 + c + t. In all of 
the equations in this paper we assume that a quark distribution function is equivalent to the corresponding 
antiparticle distribution except for u and d. In MSTW 2008, t = and it is not strictly true that s = s, 
but these are negligible effects for our calculations. Here, Gf = 1-17 x 10~ 5 GcV~ 2 is the Fermi coupling 
constant and Mn is the nucleon mass for which we use the proton mass, 0.938 GeV. The mass of the W 
boson Mw — 80.398 GeV, E v is the incident neutrino energy, and x and y are the parton momentum fraction 
(Bjorken-a;) and the inelasticity, respectively. 

Likewise, the neutral current (NC) vN cross section is given by: 
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where Mz is the Z mass. Then 
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Natural units, h = c = 1, are assumed throughout 
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FIG. 1. Calculated vN CC and NC cross sections. In this plot we compare this work (CTW), with the shaded bands 
representing the associated uncertainties due to PDFs, to those in Gandhi et al. (GQRS). Thin black lines bound 
the NC uncertainties so that they remain visible where they overlap with the CC bounds. 



and 

q o = ^ {R 2 + + * (Ll + Ll) (5) 
+ ( S + b)(Ll + Rj) + (c + t)(Ll + Rl) (6) 

with L u = 1 — 4/3-xvi/, Ld = — l + 2/3-xw, R u = — 4/3-Xw and Rd = 2/3-xw where xw — sin 2 9w = 0.226. 
For the r>N cross sections, the above equations are the same with each quark distribution function replaced 
with the corresponding antiparticlc distribution and vice versa, so that q «-> q, q° «-> q°. 

We use the parton distribution functions calculated by A.D. Martin et al. known as "MSTW 2008" [6]. 
These PDFs are the latest update to a series that began with the MRS PDFs twenty years ago, which were 
the first global next-todeading-order (NLO) PDF analysis. The MSTW 2008 set incorporates improvements 
in the precision and kinematic range of recent measurements as well as improved theoretical developments 
which make the global analysis more reliable. The publication of the MSTW 2008 set was particularly timely 
in view of the start of data taking at the Large Hadron Collider (LHC). 



B. Cross Sections 

Figs. 1 and 2 shows the results of our vN and DN cross section calculations compared to the previous 
calculations by Gandhi et al. (GQRS) [7]. These results are summarized numerically in Tables I and II. 
With regard to uncertainties, the latter paper only states that they find the uncertainties in the i>N cross 
sections to be at most a factor of 2 ±1 . Recently Cooper-Sarkar and Sarkar (CSS) [8] also published CC vN 
cross sections for energies in the range 100 < s < 10 12 GeV 2 where yfs is the COM energy of the interaction. 
There has also been a recent investigation into the dependence on the number of active heavy quarks in 
[9]. We find good agreement with the central values of both the GQRS and CSS calculations within our 
uncertainties. We note that the explicit evolution of the MSTW 2008 PDFs only takes place down to the 
lowest value of the grid points, i.e. x = 10~ 6 . Below this the values for the central set and each eigenvector 
for the error sets are extrapolated linearly in ln(l/a;). Within the region of the grids the accuracy of the 
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FIG. 2. Calculated vN CC and NC cross sections. In this plot we compare this work (CTW), with the shaded bands 
representing the associated uncertainties due to PDFs, to those in GQRS. 



NLO and NNLO evolution has been checked to small fractions of a percent, see Section 3 of [6]. 2 

Fig. 3 compares the uncertainties on the cross section calculations for the range of energies being considered. 
The uncertainties on our calculations are dramatically different from those reported by CSS for E v > 
10 8 GeV. The difference is due to a different parametrization of the gluon parton distribution g{x). The 
CSS fit to the HERA data allows a very good fit with the gluon distribution having an x dependence 
of the form g (x) oc x s . However, MSTW 2008 finds that a sum of two terms with different powers 
xg (x) oc A\x Sl + A 2 x S2 gives a better fit to the global data set. This is partially due to the global 
fit requiring a slightly larger value of the strong coupling as and consequently less gluon to drive small-x 
structure function evolution. It is also found that Tevatron jet data prefer a larger high-x gluon distribution, 
hence allowing a smaller gluon distribution at small-x from the momentum sum rule. (It is shown in [11] that 
fitting the newer combined HERA data in [12] results in no very significant change to the MSTW PDFs.) 

As well as producing the best fit, the greater flexibility brought about by including two terms in the 
parametrization results in more rapid expansion of the allowed range in g (x) at low-x beyond the reach of 
current experiments, as illustrated in Fig. 16 of [6]. With the parametrization used by CSS, the uncertainty 
can only grow as a function of ln(l/a;) in this region as discussed in Section 6.5 of [6]. We notice that 
the point at which our uncertainty starts to exceed that of CSS to a significant extent is indeed when the 
constraint due to HERA data is starting to disappear, i.e. when the dominant x values contributing to the 
cross section are x = 10 ~ 5 or lower. At very high neutrino energy when the x values probed are typically 
well below x = 10~ 5 our uncertainty on the cross section becomes very large. Hence, we conclude that 
a good measurement of the cross section within this range will provide the first direct constraint on the 
extremely small-x PDFs, and can give us the first true indication of their central value as well as reducing 
the uncertainty significantly. 

At the lower end of the uncertainty bands the cross sections do contain contributions from PDFs that 
have become negative. It is difficult to know whether this is really a problem. At low orders in as and 
leading twist perturbation theory this could lead to negative cross sections. However, we have the possibility 
of both large ln(l/x) perturbative corrections and higher- twist nonlinear effects in this regime which could 
alter this conclusion. Hence, our uncertainty at the very highest energies can be thought of as acknowledging 
the possibility of theoretical corrections in this regime. 



2 A numerical inaccuracy in the LO evolution at very small x, albeit very much less than the uncertainty, has been pointed 
out in [10]. This is unique to LO due to the extreme singular behaviour of the small-x gluon in this case. It will be corrected 
in future sets. 
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FIG. 3. Uncertainties on the calculated cross sections due to PDFs. We compare this work (CTW) to CSS. 



The recent comparative study of high energy neutrino cross sections in [13] illustrates the variation between 
models at very high energies due to theoretical assumptions. (Note the variation at lower energies between 
models in Figs. 1 and 2 of this paper is likely due to the omission of terms required at higher x [14]). 



1. Parametrizations 

For ease of use in Monte Carlo simulations, we have parametrized the cross sections in the energy range 
4 < e < 12, where e = \og 10 (E u /GcV), with an expression of the following form: 

log 10 [a(e)/cm 2 ] = d + C 2 • In (e - C ) 

+C 3 • In 2 (e - Co) (7) 
+C 4 /ln( £ -C ). 

Table III shows the values of the constants for each of uN and PN interactions, CC and NC. The 
parametrized cross sections are within approximately 1% (2%) of the calculated cross sections in the stated 
energy range for vN (PN). In Table IV we show the same constants for parametrizing the upper and lower 
bounds on the cross sections due to the uncertainties derived in this paper. For the upper bounds, the 
parametrizations are good to approximately 5% (10%) for describing our vN (PN) calculations. For the 
lower bounds, the parametrizations are at most approximately 8% from our vN (PN) calculations until 
10 u.5 < e v < 10 12 GeV where they deviate by nearly 20%. 

Note that the highest power of log 10 E v required to describe the cross section is quadratic, the same as the 
quadratic dependence of the Froissart bound [15] on hadron-hadron cross sections. This shows that although 
in principle the PDFs and cross sections grow quicker than any power of log 10 E v as E v — >• oo without some 
non-linear evolution effects slowing the evolution at very small x, in practice this has not clearly manifested 
itself in the region of energy we consider. The desire to have a parametrization for structure functions 
manifestly consistent with the Froissart bound at all energies has led to the results in [16, 17], which gives 
rather lower predictions than our central values. However, even the upper band of our uncertainty is not 
generating behaviour obviously stronger than (log 10 E v ) 2 for E v < 10 12 GeV. 

Finally, the fraction of NC events is parametrized by: 

^ = D 1 +D 2 -\n{e-D ) (8) 

"NC + &CC 



with D = 1.76, D 1 = 0.252162 and D 2 = 0.0256. 
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TABLE I. Cross Sections for z/N. 



E v (GeV) 


tree (cm 


- 2 ) 


onc (cm 


- 2 ) 


otot (cm 


- 2 ) 


1 xlO 4 


0.48x10 


-34 


0.16x10 


-34 


0.63x10 


-34 


2.5 xlO 4 


0.93x10 


-34 


0.32x10 


-34 


0.12x10 


-33 


6 xlO 4 


0.16x10 


-33 


0.57x10 


-34 


0.22x10 


-33 


1 xlO 5 


0.22x10 


-33 


0.78x10 


-34 


0.3x10" 


33 


2.5 xlO 5 


0.36x10 


-33 


0.13x10 


-33 


0.49x10 


-33 


6 xlO 5 


0.56x10 


-33 


0.21x10 


-33 


0.77x10 


-33 


1 xlO 6 


0.72x10 


-33 


0.27x10 


-33 


0.98x10 


-33 


2.5 xlO 6 


0.11x10 


-32 


0.41x10 


-33 


0.15x10 


-32 


6 xlO 6 


0.16x10 


-32 


0.61x10 


-33 


0.22x10 


-32 


1 xlO 7 


0.2xl0" : 


32 


0.76x10 


-33 


0.27x10 


-32 


2.5 xlO 7 


0.29x10 


-32 


0.11x10 


-32 


0.4x10" 


32 


6 xlO 7 


0.4xl0" : 


32 


0.16x10 


-32 


0.56x10 


-32 


1 xlO 8 


0.48x10" 


-32 


0.19x10 


-32 


0.67x10 


-32 


2.5 xlO 8 


0.67x10 


-32 


0.27x10 


-32 


0.94x10 


-32 


6 xlO 8 


0.91x10 


-32 


0.36x10 


-32 


0.13x10 


-31 


1 xlO 9 


0.11x10 


-31 


0.43x10 


-32 


0.15x10 


-31 


2.5 xlO 9 


0.14x10 


-31 


0.58x10" 


-32 


0.2x10" 


31 


6 xlO 9 


0.19x10 


-31 


0.77x10" 


-32 


0.27x10 


-31 


1 xlO 10 


0.22x10 


-31 


0.9xl0~ : 


32 


0.31x10 


-31 


2.5 xlO 10 


0.29x10 


-31 


0.12x10 


-31 


0.41x10 


-31 


6 xlO 10 


0.37x10 


-31 


0.15x10 


-31 


0.53x10 


-31 


1 xlO 11 


0.43x10 


-31 


0.18x10 


-31 


0.61x10 


-31 


2.5 xlO 11 


0.56x10 


-31 


0.23x10 


-31 


0.8x10" 


31 


6 xlO 11 


0.72x10 


-31 


0.3xl0" : 


51 


0.1x10" 


30 


1 xlO 12 


0.83x10 


-31 


0.35x10 


-31 


0.12x10 


-30 



C. Differential Cross Sections 



1. Inelasticity 

For ease of use in Monte Carlo programs, we describe here a procedure for choosing inelasticities that 
follow the proper energy-dependent distributions for energies in the range 10 7 < E v < 10 12 GeV. We use the 
Inverse Transform Method described in [18], which requires finding a function that describes dajdy which 
has an integral that is invertible. 

Due to the steepness of the differential cross section at low values of y, we divide up the parametrization 
into two regions in y: 

d* = {Y(Co,C u C 2 ) 0<y<10" 3 
dy \y'(CS,Ci) 10- 3 < y < 1 

with Y and Y' taking the following form: 

y(c ,c 1 ,c a )= g (io) 



(9) 



Y'(C>,C[) 



C' 



(ii) 
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TABLE II. Cross Sections for PN. 



E D (GeV) 


ace (cm 


- 2 ) 


onc (cm 


- 2 ) 


otot (cm 


- 2 ) 


1 xlO 4 


0.29x10 


-34 


0.11x10 


-34 


0.4x10" 


34 


2.5 xlO 4 


0.63x10 


-34 


0.24x10 


-34 


0.87x10 


-34 


6 xlO 4 


0.12x10 


-33 


0.47x10 


-34 


0.17x10 


-33 


1 xlO 5 


0.17x10 


-33 


0.67x10 


-34 


0.24x10 


-33 


2.5 xlO 5 


0.3x10" 


33 


0.12x10 


-33 


0.42x10 


-33 


6 xlO 5 


0.49x10 


-33 


0.2x10" 


33 


0.68x10 


-33 


1 xlO 6 


0.63x10 


-33 


0.26x10 


-33 


0.89x10 


-33 


2.5 xlO 6 


0.98x10 


-33 


0.4x10" 


33 


0.14x10 


-32 


6 xlO 6 


0.15x10 


-32 


0.6x10" 


33 


0.21x10 


-32 


1 xlO 7 


0.18x10 


-32 


0.76x10 


-33 


0.26x10 


-32 


2.5 xlO 7 


0.26x10 


-32 


0.11x10 


-32 


0.37x10 


-32 


6 xlO 7 


0.37x10 


-32 


0.16x10 


-32 


0.52x10 


-32 


1 xlO 8 


0.45x10 


-32 


0.19x10 


-32 


0.64x10 


-32 


2.5 xlO 8 


0.62x10 


-32 


0.27x10 


-32 


0.88x10 


-32 


6 xlO 8 


0.84x10 


-32 


0.36x10 


-32 


0.12x10 


-31 


1 xlO 9 


0.99x10 


-32 


0.43x10 


-32 


0.14x10 


-31 


2.5 xlO 9 


0.13x10 


-31 


0.58x10 


-32 


0.19x10 


-31 


6 xlO 9 


0.17x10 


-31 


0.77x10 


-32 


0.25x10 


-31 


1 xlO 10 


0.2x10" 


51 


0.9x10" 


12 


0.29x10 


-31 


2.5 xlO 10 


0.27x10 


-31 


0.12x10 


-31 


0.39x10 


-31 


6 xlO 10 


0.35x10 


-31 


0.15x10 


-31 


0.5x10" 


31 


1 xlO 11 


0.4x10" 


11 


0.18x10 


-31 


0.58x10 


-31 


2.5 xlO 11 


0.52x10 


-31 


0.23x10 


-31 


0.75x10 


-31 


6 xlO 11 


0.66x10 


-31 


0.3x10" 


11 


0.96x10 


-31 


1 xlO 12 


0.77x10 


-31 


0.35x10 


-31 


0.11x10 


-30 



TABLE III. Coefficients for parametrizing the cross sections according to Equation 7. 







Co 


Ci 


c 2 


c 3 


c 4 


V 


NC 


-1.826 


-17.31 


-6.448 


1.431 


-18.61 


V 


CC 






-6.406 




-17.91 


V 


NC 


-1.033 


-15.95 


-7.296 


1.569 


-18.30 


V 


CC 






-7.247 




-17.72 



For the low y region, the normalized integral of the distribution at y is: 

IvL Y (Co, Ci, C 2 ) dy {yo _ cr^-i/C+i) _ (ymin _ c^-V^+D 



I(Vo) = 



jy^Y(C ,C u C 2 )dy 
For the high y region, it is: 



Cl){ -i/c 2+ i) _ {yin . n _ Cl){ -i/c 2+ i)- 



i(yo) = 



J^ ia Y'(C' 0j C[)dy 
r—Y'iC^C'Jdy 



"in 


( Vo- 








In 











(12) 



(13) 



Notice that Equations 12 and 13 no longer contain the normalization factors Co and Cq. For the \ow-y 
region, y min = and j/ max = 10" 3 , while for the high-y region, j/ min = 10" 3 and y max = 1. 
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TABLE IV. Coefficients for parametrizing the uncertainty bounds on the cross sections according to Equation 7. 







C 


Ci 


c 2 


C z 


c 4 


upper 


V 


NC 


-1.456 


32.23 


-32.32 


5.881 


-49.41 




V 


cc 




33.47 


-33.02 


6.026 




V 


NC 


-2.945 


143.2 


-76.70 


11.75 


-142.8 




V 


CC 




144.5 


-77.44 


11.90 




lower 


V 


NC 


-15.35 


16.16 


37.71 


-8.801 


-253.1 




V 


CC 




13.86 


39.84 


-9.205 




V 


NC 


-13.08 


15.17 


31.19 


-7.757 


-216.1 




V 


CC 




12.48 


33.52 


-8.191 





Since Equations 12 and 13 each represent a cumulative distribution function that is invertible, we can use 
the Inverse Transform Method to select values of yo that follow the distributions in Equation 10 and 11 in 
each region. By setting /(yo) of Equations 12 and 13 to a random number R between and 1, we can then 
solve for our choice of yo and obtain for the low y region: 



yo = C\ + [R(y max - Ci)(-i/C2+i) + (i _ R ) {lJmin _ 
and for the high y region: 

(j/max — 



C* 2 /(C 2 -1) 



/yo = 



(j/min " C[) R - 



(14) 



(15) 



The parameter C\ itself depends on e, and for both regions of y, the energy dependent parameter takes 
the form: 



C 1 =A a -A 1 cxp[-(e-A 2 )/A 3 }. (16) 

In Equation 16, all parameters are primed when describing the high-y region. The numerical values of the 
parameters in Equation 16, summarized in Table V, were obtained from fits of the parametrizations to the 
theoretical calculations. 

In the low region, the parameter C 2 is also energy dependent: 

C 2 = B + B 1 ■ e (17) 

where Bo = 2.55 and B\ = —0.0949 for all interaction types. 

The fraction of the cross section occupying the low-y region is given by: 

f(s) = F -sm[F 1 -(s-F 2 )} (18) 

with ^0 = 0.128, Fi = -0.197 and F 2 = 21.8 for all interaction types. 

To summarize, for an interaction of a given type (vN or vN, CC or NC) at an energy e, one can find an 
inelasticity y chosen from the appropriate distribution through the following steps: 

• Choose a random number R\ between and 1 and if R\ < /(e) (see Equation 18), then the event lies 
in the low-?/ region. Otherwise, it is in the high-y region. 

• Obtain the value of C[, or C\ and C2, depending on the y region and the event type, using Equations 16, 
Equation 17 and Table V. 

• Choose a new random number R 2 and insert R = R 2 along with the parameters obtained in the 
previous step into Equation 14 or Equation 15 to obtain y - 
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Inelasticity y 

FIG. 4. Comparison of the uN CC inelasticity distributions from our theoretical calculation (black line) and the 
distribution obtained from the parametrization procedure described in Section II C 1 (gray histogram). 
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FIG. 5. Same as Fig. 4, but for \aw-y region. 



Figs. 4 and 5 show the calculated y distributions in each region compared with the event distributions 
generated from this procedure for vN, CC events at e = 12. 

For all event types and energies in the range 7 < e < 12, this procedure will give y distributions whose mean 
value differs from the theoretical calculation by at most 3.7% in the energy range E v = 10 7 — 10 12 GeV. Recall 
that the inelasticity is the fraction of neutrino energy carried away by the hadronic shower and therefore 
these uncertainties on the mean hadronic energy imply less than 1% uncertainties on the mean energy of the 
final state lepton. The RMS of the distributions are 0.05-0.07 in the same energy range and the difference 
in RMS values between the model and the calculation is no more than 2.6% for all interaction types except 
the vN CC events, whose uncertainties on the RMS values do not exceed 8.3%. In [19], the authors give 
a parametrization of the inelasticity distributions in the neutrino energy range 50 GeV< E v < 10 12 GeV 
which are within 15% agreement with calculations using the CTEQ6 parton distribution functions. 
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TABLE V. Coefficients that go into calculating Ci and C[ in Equation 16 for parametrizing the inelasticity distribu- 
tions. 





A 


low y 
Ai 


A 2 


A 3 




0.0 


0.0941 


4.72 


0.456 




A' 


high y 
A[ 


A' 2 


A' :i 


vN CC 
vN CC 


-0.0026 
-0.008 


0.085 
0.26 


4.1 
3.0 


1.7 


vN NC 
vN NC 


-0.005 


0.23 


3.0 


1.7 




log x 

a 10 

FIG. 6. Normalized x distribution for neutrino interactions for each decade in energy from E v = 10 4 to 10 12 GeV. 
As neutrino energy increases, the mean value of log 10 x decreases. 



Sample code for generating energy dependent inelasticity distributions for all interaction types according 
to the prescription laid out in this paper can be found at: 

http : //www. physics . ohio-state . edu/^connolly/ crosssections/y .html. 



2. Bjorken-x 

The fraction of momentum carried by a parton within a nucleon is called the Bjorkcn-x. UHE neutrino 
cross sections include contributions from PDFs in the region of x that is lower than the region above 
10~ 4 — 10~ 5 accessible in the perturbative regime by HERA experiments [12]. Fig. 6 shows the distributions 
of (1/ctot) da / d(\og 10 x) for energies between 10 4 and 10 12 GeV. As one can see the cross section starts to 
become very sensitive to the x-range below the extent of the HERA data constraint at E v ~ 10 9 GeV. 
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FIG. 7. Correlation factor cos <p between the charged current uN cross sections at the energy on the abscissa and the 
energy on the ordinate. 



D. Correlations 



Since we quote uncertainties on the vN cross sections continuously across the energy range, we also include 
the correlations between the uncertainties at different energies using the prescription laid out in [20]. For 
completeness we briefly summarize the procedure here. 

Each cross section value calculated in this paper is the sum of contributions from TV different orthogonal 
eigenvectors which are the result of diagonalizing the parameters of the PDFs. Considering two different 
cross sections X and Y, the correlation between their uncertainties is denoted cos ip and is given by: 



COSp = 



AX-AY 
AXAY 
1 

4AXAY 



(19) 



where 



AX 



AX 



1 



JV 



E(*< (+) -** (_) 



i=l 



(20) 



Here, x[ +) and x\ ] are the upper and lower bounds on the contribution to X from the i th eigenvector. 
For a maximum correlation, cos tp — 1, for an anticorrelation, coscp = — 1 and for two quantities that are 
uncorrelated, cos ip = 0. 

In Fig. 7, we plot cos ip for all energy pairs in the range 10 4 — 10 12 GeV for charged current vN cross 
sections. The analogous plots for vN and neutral current cross sections look similar. There are very strong 
correlations among cross sections at energies above 10 9 GeV, and then again below 10 7 GeV, with little 
correlation between energies in different regions. There are three reasons for this, all associated with the 
fact that in the lower energy region the cross section has a high proportion of its contribution from x <~ 0.01 
or above, while in the higher energy region most of the contribution is from x values lower than this. First, 
from sum rules in the PDFs, there is a crossing point where changes in PDFs become anticorrelated, i.e. for 
any change the PDFs, they tend to increase below this x and decrease above this x or vice versa. For high 
energy scales this is at x w 0.01. Second, this also happens to be the x where there is a very large amount 
of accurate HERA data also tending to fix the PDFs. Thirdly, x — 0.01 is a transition point at which for 
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higher x the dominant contributions are from valence quarks but for lower x the sea quarks dominate, which 
are gluon driven. 

III. CROSS SECTION CONSTRAINTS 
A. Motivation 

The upper (lower) bounds on the Standard Model vN cross sections differ from the central values by more 
than approximately 20% (15%) for neutrino energies above 10 9 GeV. Above 10 10 GeV, the uncertainties are 
more than approximately 50% (30%). We would like to constrain the cross sections at the highest energies 
because in that region neutrino experiments could be sensitive to new physics scenarios. Fortunately, there 
is a near guaranteed flux of neutrinos in the UHE energy regime from GZK interactions, but this so-called 
"cosmogenic" neutrino flux has large theoretical uncertainties associated with it. 

Here we outline a technique [21] to constrain the UHE vN cross sections that is independent of the 
incident neutrino flux through the measured zenith angle distributions with a subterranean detector such 
as IceCube, ARA or ARIANNA [22-24]. The latter two experiments are currently in the first stages of 
deploying prototype detectors and will be focused on the UHE regime. Here, we focus on the ARA detector 
as an example, but our results are general to any subterranean neutrino detector with similar capabilities in 
energy measurement and reconstruction. 

The sensitivity of neutrino detectors to uN cross sections due to earth absorption has been addressed 
elsewhere in varying degrees [25-33]. This is the only study that uses the full zenith angle distribution in 
the UHE regime to make hard predictions for the expected constraints on models with extra dimensions. 

We note a few assumptions made for this study. First, we assume the incident flux is entirely neutrinos with 
no anti-neutrino component. Above 10 8 GeV, the cross sections for neutrinos and anti-neutrinos differ by 
no more than about 6%. Second, we assume that all of the interactions occur at the same depth, d = 250 m, 
and that it is precisely known. In an actual data analysis one would modify the dn/d cos 9 Z distributions 
so that d is the measured depth, just as we use the measured energy for each event. The three-dimensional 
vertex resolution of the ARA detector is expected to be of order 10's of meters. One can show that a depth 
uncertainty given by 5d changes dP/dcos9 z by a fraction of order 5d/L. This only approaches of order 
10% at the highest cross sections probed in this paper, approximately 10~ 28 - 5 cm 2 . We also assume that 
the detector efficiency is flat in zenith angle 9 Z and that the neutrino flux is isotropic. An analysis would 
fold in the measured detector efficiency, and neutrinos from any discovered sources could be removed from 
the sample. Finally, we assume that the energy resolution A log 10 E v = 0.4, which is dominated by the 
uncertainty in inelasticity, is appropriate for both neutrino and black hole interactions. 

B. ARA 

The Askaryan Radio Array (ARA) is a nascent neutrino detector near the South Pole. It is an array of radio 
antennas deployed deep in the ice, designed to measure the radio Cerenkov pulse from UHE neutrinos [34- 
37]. If expanded to become a precision measurement, observatory class array of 300 to 1000 km 2 area, ARA 
will be capable of measuring hundreds of cosmogenic neutrinos per year. The first ARA testbed station was 
deployed in the 20102011 austral summer and the first ARA stations will be deployed in 2011-2012. 

For the purpose of generating mock ARA data for this study, we use energy- dependent relative effective 
areas (the overall scale is not used) derived by inverting the projected flux limits for a 37-station array in [23] . 
The stations are arranged on a triangular grid with the array forming a hexagon. Each station is made of 
three "strings" deployed vertically in the ice, each holding two pairs of vertically and horizontally polarized 
antennas that sit at 200 m depth, for a total of 12 antennas per station. The trigger requires 5 out of 12 
stations measure a pulse that exceeds 3.5 times the expected noise level. We assume an energy resolution 
of Alog 10 _E l/ = 0.4. The resolution on the neutrino 9 Z is expected to be approximately 2° for the events 
(approximately 80%) detected by one station only [38]. 

Our projected constraints will not depend on the exact size of the ARA detector. This is because the 
expected limits will be quoted for a specific number of neutrino events measured in the detector, whether 
that came from a weak flux measured with a large detector or a strong flux measured with a smaller detector. 
Nor do the projected constraints depend strongly on the energy threshold of the experiment, as long as there 
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is at least a crude energy measurement capability. We will quantify these statements at the conclusion of 
Section III E. In addition, our conclusions are not specific to the ARA experiment, and would be qualitatively 
similar for any subterranean detector with sensitivity to UHE neutrinos. 



C. A Flux-Independent Technique for Measuring the UHE Neutrino-Nucleon Cross Section 

Consider a neutrino (see Figure 8) that interacts in a subterranean detector such as ARA at depth d and 
zenith angle 6 Z , the angle from vertical of the direction of origin of each incident neutrino. The neutrino 
travels a distance D through the earth of radius R before reaching its interaction point. From Figure 8 and 
using the law of sines, 



sin0i sin(7r — 9 Z ) 
R-d ~ R ' 

which gives 



(21) 



R 

Then since 02 = n/2 — 0i, 



R — d . 

sin 0i = — - — smO z . (22) 



sin 4>2 — cos 0i = \J 1 — sin 2 0i = y 1 — — ^ sin 2 9 Z (23) 

Then the distance traveled through the earth by the neutrino is: 

D = (D - x) + x = #sin0 2 + (R - d) sin (e z - |) (24) 
Inserting Eq. 23 into Eq. 24 we find: 



D = R ^ 1 ~(^1T~) sin2 ^ + (^- rf ) sin (^-f) (25) 
Then, taking d « R and replacing sin# 2 with \/l — cos 2 9 Z , we find 



D = yj (R 2 - 2Rd) cos 2 9 Z + 2Rd -(R-d) cos 9 Z . (26) 

If the detection efficiency is uniform in 9 Z , then for a neutrino with energy E v , the probability distribution 
in 9 Z for a detected interaction is given by [39]: 



dP m =A , cxp ,_D\ =A (_^_ 2Rd) cos 2 + 2Rd -(R-d) cos9 z 



dcos9 z y "V L(E V )J 1 \ L(E V ) 

where R is the radius of the earth, A is a constant that sets the total probability to unity, and L (E v ) is the 
interaction length for a neutrino of energy E v along its path through the earth, given by: 

°(^ v ) (p) 0z 

where Mn is the nucleon mass, <j(E v ) is the vN cross section at energy E v and (p) s is the mean density 
averaged over the distance travelled by the neutrino at the given 9 Z . The expected 9 Z distribution from 
a sample of n measured neutrino interactions with measured energies E z v will be the sum of n different 
(dP/dcos9 z )i, so that the resulting expected 9 Z distribution is then given by: 



dn A ( V '( R2 - 2Rd ) cos2 9 z +2Rd- (R-d) cos 9 Z , 

— y, Ai ■ cxp / ~ \ • (29) 

- 1 L[El) 



dcos9 zcxp 
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FIG. 8. Diagram showing neutrino incident on the earth used in derivation of Eq. 27. 




FIG. 9. Probability distributions in zenith angle for select neutrino energies. The width of each band is due to the 
cross section uncertainties reported in this paper. The kinks in the distributions at cos# z = —0.1 and cosf^ = —0.8 
are due to the neutrino paths reaching the earth's mantle and core, respectively. 



Note that due to the difference between measured neutrino energies E l v and true energies ££, the expected 
distribution will differ from the true one. In Fig. 9, we plot dP/dcos9 z for a few monoenergetic distributions. 
Notice the breaks in the distribution due to the neutrino trajectory intersecting the Earth's core and mantle, 
derived using a simple three-layer model of the Earth's interior. Due to this structure at lower energies, this 
technique could lead to a measurement of the Earth's density profile with a km 3 neutrino experiment such 
as IceCube that is independent of the traditional techniques used by geologists [39] . 
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log ( E / GeV ) 

FIG. 10. Cross sections for uN interactions in models with extra space-time dimensions compared with the SM vN 
cross sections. The gray band surrounding the SM cross sections are the uncertainties presented in this paper. 



D. Models with Extra Space-Time Dimensions 



There are a class of models for physics beyond the Standard Model that contain extra space-time dimen- 
sions [40-43] . These models are motivated by the need to resolve what is known as the Hierarchy Problem 
in particle physics, where the dramatically different energy scales for electroweak symmetry breaking and 
quantum gravity lead to a need for fine tuning of terms in the calculation of the Higgs mass. 

In these extra-dimensional models (EDMs), the energy scale at which gravity dominates, Mjj, is reduced 
to of order 1 TeV, just above the electroweak scale. The weakness of gravity in our 3+1 dimensional world is 
a consequence of its propagation in additional dimensions. The number of dimensions in the model beyond 
the four known space-time dimensions is denoted Np- 

Interactions at energies at or above the reduced Planck mass lead to the production of micro-black holes, 
and this additional channel causes cross sections to be enhanced. The minimum black hole mass is given by 
^sh ~ XminMn, where x min is a parameter in the model. Fig. 10 shows the predicted vN cross sections 
for a few EDMs, from [43], compared to the SM cross sections calculated in this paper. 

Tcvatron experiments CDF and DO have already set lower limits on Md in the range of approximately 
1-1.6 TeV with between 2 and 7 extra dimensions [44, 45]. A recent paper by the CMS collaboration places 
constraints on black hole production at the LHC based on 35 pb _1 of data at center-of-mass energy of 
7 TeV [46]. For x m ; n = 1, they exclude models with n up to 6 for 1.5 < Md < 3.5 TeV. They were not 
sensitive to models with a; m ; n = 3 in the range of Md considered. In addition, neutrino experiments have 
constrained UHE cross sections, but in a way that depends on a model for the neutrino flux [47]. With the 
technique described here, neutrino experiments would be setting limits on EDMs that are competitive with 
collider experiments and independent of a flux model. 

In Fig. 11, we plot the expected Z distribution for 100 events measured in ARA using cross sections from 
the SM compared to ones from a few select EDMs and a bin width of A cos 6* z = 0.1. The shape of the 
true energy spectrum of the 100 events is the product of an incident flux spectrum and an energy dependent 
effective area for the detector. Here we assume a neutrino flux from the GZK process as in [48], and the 
ARA effective areas derived from [23] as described in Section IIIB. 
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FIG. 11. Predicted cos9 z distributions for 100 neutrino events observed by ARA for the SM and the four different 
EDMs shown in Figure 10.. The bin width is taken to be Acos6> z = 0.1. 



E. Projected Constraints 

In order to assess the sensitivity of a future ARA detector to EDMs, we generate many pseudoexper- 
imcnts, and compare the resulting pseudodata distributions to the predicted ones for the signal and null 
hypothesis respectively. For an expected number of events N exp , the number of events observed in a given 
pseudoexperiment is given by n p and is Poisson distributed with mean A cxp . The data in the i th bin in 
cos Z is denoted n Pj i. The mean number of events predicted in the i th bin centered on cos6* z = co is 

rfri 

VM-Acos^ (30) 



dcosO 



where dn p /d cos 8 Z is constructed for each pseudoexperiment using Equation 29. Note that through the 
L(££)'s in Equation 29, dn p /dcos6 z depends on the model hypothesis, whether that be the SM or an EDM. 
Also note that we only use the measured energies, selected from the incident flux spectrum and then smeared 
according to detector resolution of A\ogE v = 0.4. The bin width is Acos6* z = 0.1, or A8 Z w 5.7°, which is 
greater than the expected resolution of the ARA array in 9 Z . 

We use the following ratio of Poisson probabilities to discriminate between the two hypotheses [49] : 

Q _ P P oi ss (data | EDM truth) 

W P poiss (data | SM truth) [ ' 



where 



N n 11 ' g-fiEDM,i 

P poiss (data | EDM truth) = JJ EDM ' ! (32) 

i=i ni ' 



and 



P poiss (data | SM truth) = JJ bM '' ■ . (33) 

i=l n '' 

Here, N is the number of bins, and n, is the number of events measured in the i th bin. The number of events 
expected in a bin from an extra-dimensional model is /iedm and the number expected in the Standard Model 
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FIG. 12. Comparison of likelihood ratio — 21nQ from iVp pseudoexperiments generated assuming SM and EDM 
truth. The region between — 21nQi 6 and — 21nQg4 contains 68% of the SM pseudoexperiments, with the median at 
-21nQ 50 . 



is Msm- Then, wc find 



The parameters for the EDM models are denned in Section HID. Equation 34 is evaluated separately 
for pseudodata n p generated assuming SM and EDM truths, giving a different — 21nQ distribution for 
each. According to the Neyman-Pearson lemma, this likelihood ratio is the test statistic with the most 
discriminating power [50]. 

We estimate the constraints on EDMs expected to be set by the full scale ARA detector described in 
Section IIIB depending on the observed number of events n p . For the expected limit, we first consider the 
median value — 21nQ5o of the — 21nQ distribution from SM pseudoexperiments. The subscript denotes the 
percentage of SM pseudoexperiments with lower values of — 21nQ. Then, for a given EDM model, the p 
value is the fraction of EDM pseudoexperiments with —2 In Q < —2 In Q50. Then, on average we can expect 
that a model can be excluded with percentage confidence level CL = 100 x (1 — p). 

The observed limit will differ from that expected due to fluctuations in the data, as reflected by the width 
of the —2 InQ distribution for SM pseudoexperiments. Therefore, we quote a range for the expected CLs by 
calculating p values for — 21nQi6 and — 21n<5s4, so that 68% of pseudoexperiments would give confidence 
levels in the range. Fig. 12 shows the distributions in —2 InQ for the SM and one EDM hypothesis, for 100 
events observed in ARA. We mark the points along the abscissa — 21nQi6, — 21nQ5o and — 21n<5s4- 

Fig. 13 shows the expected CLs for constraining EDMs as a function of the number of events observed with 
ARA. The black bands show the variation in the expected constraints brought about by the uncertainties 
on the SM cross sections presented here. The gray bands show the range of expected limits due to variations 
in the data; sixty-eight percent of the pseudoexperiments give constraints that lie in the gray region. The 
gray bands are centered on the curve corresponding to the central value standard model cross sections. 
Conservatively taking the SM upper bounds to be the true cross sections, for 100 events observed with ARA, 
the mean expectation is to exclude the following models: x m j n = 1, Md = 1, Ad > 2; x m ; n = 3, Md = 
1, Ad > 3; x m ; n — 1, Md = 2, TVd > 3. For = 3, Md = 2, 110 events would be needed to be 

predicted to exclude Ad = 7. 

We have checked that these projected limits are robust to changes in the neutrino spectrum as well as 
details of the sensitivity of the experiment. We have calculated the expected constraints with an E~ 3 
neutrino spectrum, and also after shifting the ARA effective areas up and down in energy by a factor of 3. 
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FIG. 13. Predicted confidence levels for excluding selected EDMs with an ARA experiment, depending on the 
number of events observed. The upper (lower) edges of the black bands in the figure show the range of expected 
limits taking the SM cross section to be their lower (upper) bounds as presented in this paper. The gray bands show 
the uncertainty due to pseudoexperiment statistics, with 68% of experiments expected to result in limits that lie 
within the gray region (using the central values for the SM cross section as the true cross sections). 



All of these changes have less of an effect on the expected constraints than the cross section uncertainties 
(the effect of which are depicted in the black bands of Fig. 13). 



IV. CONCLUSIONS 



We have presented new calculations of the vN and vN , CC and NC cross sections in the neutrino energy 
range 10 4 < E v < 10 12 GcV using the MSTW 2008 PDFs, along with their PDF uncertainties. The cross 
section values are consistent with those reported in previous publications but the uncertainties presented here 
are significantly larger. This difference is due to a two-parameter model for the gluon parton distribution 
at very small x used by MSTW 2008 arising from a fit to the global data, which leads to a wide range of 
allowed values for the gluon contribution at low- a;. For ease of use in Monte Carlo simulations, we have 
provided parametrizations of the cross sections for each event type, the fraction of each type and the cross 
section uncertainty bounds in the energy range. In addition, we have outlined a procedure for generating y 
values for a neutrino sample with E v > 10 7 GeV using the Inverse Transform Method. 
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Finally, we present a technique for constraining UHE cross sections with a next-generation subterranean 
neutrino experiment. Using ARA as an example, we have shown that with 100 events observed, neutrino 
experiments set constraints on extra-dimensional models that arc competitive with those set by collider 
experiments. 
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